% Written as if Fixed Cost NOT Subsidized !!!!!
% ====================================================================

clear all; clear global; clc; close all;

load('../Set_postestimation_use.mat')

tau_A = 0.053;
tau_B = 0.038;

tau_bA = 0.053;
tau_bB = 0.041;

Par.tau_A = tau_A;
Par.tau_B = tau_B;

Par.L_A   = 1;
Par.L_B   = 1;

Par.time_ss  = 2*Par.time_ss;
Par.time_81  = Par.time1;
Par.time_95  = 20/Par.dt;
Par.time_add = (200+Par.time_81*Par.dt)/Par.dt-Par.time_ss+1;

[Par,Step,Trans] = Trade_ctrf(Par);

horizons = [5:5:50];
tau_bA_grid = linspace(0.00,0.70,151);
trfA_grid = 0.10:-0.01:0;


rev_b   = 1; % Revenue back to consumer
retal_b = 0; % 1 if foreign retaliation

Par.rev_US = rev_b;
Par.rev_FN = rev_b;

% ============================================
%                 SOLUTIONS
% ============================================

Param = structvars(Par);
for ctr0 = 1:size(Param,1); eval(Param(ctr0,:)); end

steps      = Step{1};
Trans_innov = Step{2};
Trans_fail = Step{3};
Trans_exo  = Step{4};

%% Calibrated Economy Path
% =================================================================

Tol_rr = 1*10^-6;
save Tol_rr Tol_rr

%---------------- Path

mu_init = y_init;
Q_init  = y_init;

Par.if_long = 0;

[rr_sqn, X_sqn, V_sqn, Xe_sqn, Q_sqn, Ql_sqn, mu_share, ind_neg, ...
                            prft_fnc_sqn, wage_dom_sqn, mcuts] = ...
                                Path_ctrf(mu_init,Q_init,Par,Step,Trans);

mu_init = mu_share(:,Par.time_81);
Q_init  = Q_sqn(1:size(Q_sqn,1)/2,Par.time_81);

Par.time_81 = 1;

parpool(20)

for ctr_k = 1:length(trfA_grid)

    Par.trf_A = trfA_grid(ctr_k);  
    Par.trf_B = trfA_grid(ctr_k); 

    [Par,Step,Trans] = Trade_ctrf(Par);
    
    Param = structvars(Par);
    for ctr0 = 1:size(Param,1); eval(Param(ctr0,:)); end
    
    steps      = Step{1};
    Trans_innov = Step{2};
    Trans_fail = Step{3};
    Trans_exo  = Step{4};

    Tol_rr = 1*10^-6;
    save Tol_rr Tol_rr
    
    [rr_sqn, X_sqn, V_sqn, Xe_sqn, Q_sqn, Ql_sqn, mu_share, ind_neg, ...
                                prft_fnc_sqn, wage_dom_sqn, mcuts] = ...
                                    Path_ctrf(mu_init,Q_init,Par,Step,Trans);
                                
    %---------------- Variables
    
    % Profits
    % --------------------------------------------------------
    
    prft_fnc_A = prft_fnc_sqn(1:size(prft_fnc_sqn,1)/2,:);
    prft_fnc_B = prft_fnc_sqn(size(prft_fnc_sqn,1)/2+1:end,:);
    
    wage_dom_A = wage_dom_sqn(1:size(wage_dom_sqn,1)/2,:);
    wage_dom_B = wage_dom_sqn(size(wage_dom_sqn,1)/2+1:end,:);
    
    % Qualities
    % --------------------------------------------------------
    
    QA = Q_sqn(1:size(Q_sqn,1)/2,:);
    QB = Q_sqn(size(Q_sqn,1)/2+1:end,:);
    
    QA_long = Ql_sqn(1:size(Q_sqn,1)/2,:);
    QB_long = Ql_sqn(size(Q_sqn,1)/2+1:end,:);
    
    QwA = sum(wage_dom_A.*QA+flip(1-wage_dom_B).*QA*(1+kapa)/((1+kapa)*(1+trf_B))^(1/bet));
    QwB = sum(wage_dom_B.*QB+flip(1-wage_dom_A).*QB*(1+kapa)/((1+kapa)*(1+trf_A))^(1/bet));
    
    qq_sqn = [sum(QA)./QwA; sum(QB)./QwB];
                    
    wq_sqn = [(1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn(1,:).^(-Par.bet);
                (1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn(2,:).^(-Par.bet)];  
    
    xk_sqn = [(Par.pi_frn/Par.bet)^(1/(1-Par.bet))./wq_sqn(1,:).^(1/Par.bet)*Par.L_B;
             (Par.pi_frn_toUS/Par.bet)^(1/(1-Par.bet))./wq_sqn(2,:).^(1/Par.bet)*Par.L_A];
    
    % Shares
    % --------------------------------------------------------
    
    MU      = mu_share(:,1:time_ss+1);
    MU_long = mu_share;
    
    McutX  = mcuts(1,1:time_ss);
    McutM  = mcuts(2,1:time_ss);
    
    % R&D
    % --------------------------------------------------------
    
    X_iA = X_sqn(1:size(X_sqn,1)/2,:);
    X_iB = X_sqn(size(X_sqn,1)/2+1:end,:);
    
    X_eA = Xe_sqn(1:size(Xe_sqn,1)/2,:);
    X_eB = Xe_sqn(size(Xe_sqn,1)/2+1:end,:);
    
    % 2) Income 
    % --------------------------------------------------------
    
    II.IA_1_pr  = sum(prft_fnc_A.*QA).*wq_sqn(1,:).^(1-1/bet);
    II.IA_2_rdi = sum(alfa_A/gama_A*X_iA.^gama_A.*QA);
    II.IA_3_rde = sum(alfa_eA/gama_eA*X_eA.^gama_eA.*QA);
    II.IA_4_wd  = pi_dom/(1-bet)*L_A*wq_sqn(1,:).^(1-1/bet).*sum(wage_dom_A.*QA);
    II.IA_5_wf  = pi_frn_toUS/(1-bet)*L_A*wq_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_dom_A).*QB);
    II.IA_6_wig = wq_sqn(1,:).*sum(QA);
    II.IA_7_trf = rev_US*trf_A*((1+trf_A)*(1+kapa)*eta/(1-bet))*wq_sqn(2,:).*xk_sqn(2,:).*sum(flip(1-wage_dom_A).*QB);
    
    II.IB_1_pr  = sum(prft_fnc_B.*QB).*wq_sqn(2,:).^(1-1/bet);u
    II.IB_2_rdi = sum(alfa_B/gama_B*X_iB.^gama_B.*QB);
    II.IB_3_rde = sum(alfa_eB/gama_eB*X_eB.^gama_eB.*QB);
    II.IB_4_wd  = pi_dom/(1-bet)*L_B*wq_sqn(2,:).^(1-1/bet).*sum(wage_dom_B.*QB);
    II.IB_5_wf  = pi_frn/(1-bet)*L_B*wq_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_dom_B).*QA);
    II.IB_6_wig = wq_sqn(2,:).*sum(QB);
    II.IB_7_trf = rev_FN*trf_B*((1+trf_B)*(1+kapa)*eta/(1-bet))*wq_sqn(1,:).*xk_sqn(1,:).*sum(flip(1-wage_dom_B).*QA);
    
    IA = sum(prft_fnc_A.*QA).*wq_sqn(1,:).^(1-1/bet)-...
        sum(alfa_A/gama_A*X_iA.^gama_A.*QA)-...
        sum(alfa_eA/gama_eA*X_eA.^gama_eA.*QA)+...
        pi_dom/(1-bet)*L_A*wq_sqn(1,:).^(1-1/bet).*sum(wage_dom_A.*QA)+...
        pi_frn_toUS/(1-bet)*L_A*wq_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_dom_A).*QB)+...
        wq_sqn(1,:).*sum(QA)+...
        rev_US*trf_A*((1+trf_A)*(1+kapa)*eta/(1-bet))*wq_sqn(2,:).*xk_sqn(2,:).*sum(flip(1-wage_dom_A).*QB);  
    IB = sum(prft_fnc_B.*QB).*wq_sqn(2,:).^(1-1/bet)-...
        sum(alfa_B/gama_B*X_iB.^gama_B.*QB)-...
        sum(alfa_eB/gama_eB*X_eB.^gama_eB.*QB)+...
        pi_dom/(1-bet)*L_B*wq_sqn(2,:).^(1-1/bet).*sum(wage_dom_B.*QB)+...
        pi_frn/(1-bet)*L_B*wq_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_dom_B).*QA)+...
        wq_sqn(2,:).*sum(QB)+...
        rev_FN*trf_B*((1+trf_B)*(1+kapa)*eta/(1-bet))*wq_sqn(1,:).*xk_sqn(1,:).*sum(flip(1-wage_dom_B).*QA); 
    
    g_IA = mean((IA(2:time_81+1)-IA(1:time_81))./IA(1:time_81));
    g_IB = mean((IB(2:time_81+1)-IB(1:time_81))./IB(1:time_81)); 
    
    gI_A = (IA(2:end)./IA(1:end-1))-1;
    gI_B = (IB(2:end)./IB(1:end-1))-1;
    
    % Welfare
    % --------------------------------------------------------
    
    if psy == 1
        
        for ctr_t = 1:length(horizons)
        
            Wlf_o(ctr_t,1) = sum(exp(-ro*(dt:dt:horizons(ctr_t))).*log(IA(time_81+1:time_81+horizons(ctr_t)/dt)));
    
            Wlf_oB(ctr_t,1) = sum(exp(-ro*(dt:dt:horizons(ctr_t))).*log(IB(time_81+1:time_81+horizons(ctr_t)/dt)));
            
        end  
        
    else
        
        for ctr_t = 1:length(horizons)
            
            Wlf_o(ctr_t,1) = (1-psy)^-1*[sum(exp(-ro*(dt:dt:horizons(ctr_t))).*(IA(time_81+1:time_81+horizons(ctr_t)/dt).^(1-psy)))];
    
            Wlf_oB(ctr_t,1) = (1-psy)^-1*[sum(exp(-ro*(dt:dt:horizons(ctr_t))).*(IB(time_81+1:time_81+horizons(ctr_t)/dt).^(1-psy)))];
            
        end  
       
    end
      
    
    % Policy change
    % =================================================================
    %
    
    % ------------- Parfor intialization ------------- 
    
    mub_init = mu_share(:,time_81+1);
    Qb_init  = QA(:,time_81+1);
    
    mu_pre = mu_share(:,1:time_81);
    rr_pre = rr_sqn(:,1:time_81);
    IA_pre = IA(:,1:time_81);
    IB_pre = IB(:,1:time_81);
    QA_pre = QA(:,1:time_81);
    QB_pre = QB(:,1:time_81);
    QwA_pre = QwA(:,1:time_81);
    QwB_pre = QwB(:,1:time_81);
    McutM_pre = McutM(1,1:time_81);
    McutX_pre = McutX(1,1:time_81);
    
    QA_post = QA(:,time_81+1:end);
    QB_post = QB(:,time_81+1:end);
    
    II_orig = II;
    trfA_call = trfA_grid(ctr_k);
    
    Tol_rr = 5*10^-5;
    save Tol_rr Tol_rr

    parfor ctr_s = 1:length(tau_bA_grid)
        
        Run = [ctr_s,ctr_k]

        Par_b = Par;
        Par_b.num_r = Par.num_r/2;

        Par_b.taup_A = 0;
        Par_b.taup_B = 0;

        Par_b.tau_A = tau_bA_grid(ctr_s);
        Par_b.tau_B = tau_bB;

        Par_b.tau_eA = 0.0;
        Par_b.tau_eB = 0.0;

        [Par_b,Step_b,Trans_b] = Trade_ctrf(Par_b);

        % ============================================
        %                 SOLUTIONS
        % ============================================

        pi_domb = Par_b.pi_dom;
        pi_frnb = Par_b.pi_frn;
        pi_frnb_toUS = Par_b.pi_frn_toUS;

        %---------------- Path

        [rrb_sqn, Xb_sqn, Vb_sqn, Xeb_sqn, Qb_sqn, Qlb_sqn, mub_share, ind_neg_b, ...
                                    prft_fncb_sqn, wage_domb_sqn, mcutsb] = ...
                                    Path_ctrf(mub_init,Qb_init,Par_b,Step_b,Trans_b);

        rr_b = [rr_pre rrb_sqn(:,1:end-time_81)];

        %---------------- Variables

        % Profits
        % --------------------------------------------------------

        prft_fncb_A = prft_fncb_sqn(1:size(prft_fncb_sqn,1)/2,:);
        prft_fncb_B = prft_fncb_sqn(size(prft_fncb_sqn,1)/2+1:end,:);

        wage_domb_A = wage_domb_sqn(1:size(wage_domb_sqn,1)/2,:);
        wage_domb_B = wage_domb_sqn(size(wage_domb_sqn,1)/2+1:end,:);

        % Qualities
        % --------------------------------------------------------

        QA_b = Qb_sqn(1:size(Q_sqn,1)/2,:);
        QB_b = Qb_sqn(size(Q_sqn,1)/2+1:end,:);

        QAb_long = Qlb_sqn(1:size(Q_sqn,1)/2,:);
        QBb_long = Qlb_sqn(size(Q_sqn,1)/2+1:end,:);

        QwA_b = sum(wage_domb_A.*QA_b+flip(1-wage_domb_B).*QA_b*(1+kapa)/((1+kapa)*(1+Par_b.trf_B))^(1/bet));
        QwB_b = sum(wage_domb_B.*QB_b+flip(1-wage_domb_A).*QB_b*(1+kapa)/((1+kapa)*(1+Par_b.trf_A))^(1/bet));

        qqb_sqn = [sum(QA_b)./QwA_b; sum(QB_b)./QwB_b];

        wqb_sqn = [(1-Par_b.bet)*Par_b.eta^(Par_b.bet-1)*qqb_sqn(1,:).^(-Par_b.bet);
                    (1-Par_b.bet)*Par_b.eta^(Par_b.bet-1)*qqb_sqn(2,:).^(-Par_b.bet)];  

        xkb_sqn = [(Par_b.pi_frn/Par_b.bet)^(1/(1-Par_b.bet))./wqb_sqn(1,:).^(1/Par_b.bet)*Par.L_B;
                    (Par_b.pi_frn_toUS/Par_b.bet)^(1/(1-Par_b.bet))./wqb_sqn(2,:).^(1/Par_b.bet)*Par.L_A];

        QA_orgb = [QA_pre QA_b(:,1:end-time_81)];
        QB_orgb = [QB_pre QB_b(:,1:end-time_81)];

        QwA_orgb = [QwA_pre QwA_b(:,1:end-time_81)];
        QwB_orgb = [QwB_pre QwB_b(:,1:end-time_81)];

        qq_orgb  = [sum(QA_orgb)./QwA_orgb; sum(QB_orgb)./QwB_orgb];

        % Shares
        % --------------------------------------------------------

        MU_b = [mu_pre mub_share(:,1:end-time_81)];

        MCUTX(ctr_k,ctr_s).cut  = [McutX_pre mcutsb(1,1:time_ss-time_81)];
        MCUTM(ctr_k,ctr_s).cut  = [McutM_pre mcutsb(2,1:time_ss-time_81)];

        % R&D
        % --------------------------------------------------------

        Xb_iA = Xb_sqn(1:size(X_sqn,1)/2,:);
        Xb_iB = Xb_sqn(size(X_sqn,1)/2+1:end,:);

        Xb_eA = Xeb_sqn(1:size(Xe_sqn,1)/2,:);
        Xb_eB = Xeb_sqn(size(Xe_sqn,1)/2+1:end,:);

        % Income 
        % --------------------------------------------------------

        II = II_orig;

        II.IA_bb1_pr  = sum(prft_fncb_A.*QA_b).*wqb_sqn(1,:).^(1-1/bet);
        II.IA_bb2_rdi = sum(alfa_A/gama_A*Xb_iA.^gama_A.*QA_b);
        II.IA_bb3_rde = sum(alfa_eA/gama_eA*Xb_eA.^gama_eA.*QA_b);
        II.IA_bb4_wd  = pi_domb/(1-bet)*L_A*wqb_sqn(1,:).^(1-1/bet).*sum(wage_domb_A.*QA_b);
        II.IA_bb5_wf  = pi_frnb_toUS/(1-bet)*L_A*wqb_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_domb_A).*QB_b);
        II.IA_bb6_wig = wqb_sqn(1,:).*sum(QA_b);
        II.IA_bb7_trf = Par_b.rev_US*Par_b.trf_A*((1+Par_b.trf_A)*(1+kapa)*eta/(1-bet))*wqb_sqn(2,:).*xkb_sqn(2,:).*sum(flip(1-wage_domb_A).*QB_b);

        II.IA_bb2_xi = sum(alfa_A/gama_A*X_iA.^gama_A.*[QA(:,1:time_81) QA_b(:,1:end-time_81)]);
        II.IA_bb3_xe = sum(alfa_eA/gama_eA*X_eA.^gama_eA.*[QA(:,1:time_81) QA_b(:,1:end-time_81)]);

        II.IB_bb1_pr  = sum(prft_fncb_B.*QB_b).*wqb_sqn(2,:).^(1-1/bet);
        II.IB_bb2_rdi = sum(alfa_B/gama_B*Xb_iB.^gama_B.*QB_b);
        II.IB_bb3_rde = sum(alfa_eB/gama_eB*Xb_eB.^gama_eB.*QB_b);
        II.IB_bb4_wd  = pi_domb/(1-bet)*L_B*wqb_sqn(2,:).^(1-1/bet).*sum(wage_domb_B.*QB_b);
        II.IB_bb5_wf  = pi_frnb/(1-bet)*L_B*wqb_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_domb_B).*QA_b);
        II.IB_bb6_wig = wqb_sqn(2,:).*sum(QB_b);
        II.IB_bb7_trf = Par_b.rev_FN*Par_b.trf_B*((1+Par_b.trf_B)*(1+kapa)*eta/(1-bet))*wqb_sqn(1,:).*xkb_sqn(1,:).*sum(flip(1-wage_domb_B).*QA_b);

        II.IB_bb2_xi = sum(alfa_B/gama_B*X_iB.^gama_B.*[QB(:,1:time_81) QB_b(:,1:end-time_81)]);
        II.IB_bb3_xe = sum(alfa_eB/gama_eB*X_eB.^gama_eB.*[QB(:,1:time_81) QB_b(:,1:end-time_81)]);

        IA_bb = II.IA_bb1_pr-II.IA_bb2_rdi-II.IA_bb3_rde+...
                II.IA_bb4_wd+II.IA_bb5_wf+II.IA_bb6_wig+II.IA_bb7_trf;

        IB_bb = II.IB_bb1_pr-II.IB_bb2_rdi-II.IB_bb3_rde+...
                II.IB_bb4_wd+II.IB_bb5_wf+II.IB_bb6_wig+II.IB_bb7_trf;

        IA_b = [IA_pre IA_bb(:,1:end-time_81)];
        IB_b = [IB_pre IB_bb(:,1:end-time_81)];

        % Alternative for Welfare

        QA_bo = [QA_post 0*QA_b(:,end-time_81+1:end)];
        QB_bo = [QB_post 0*QB_b(:,end-time_81+1:end)];

        QwA_bo = sum(wage_domb_A.*QA_bo+flip(1-wage_domb_B).*QA_bo*(1+kapa)/((1+kapa)*(1+Par_b.trf_B))^(1/bet));
        QwB_bo = sum(wage_domb_B.*QB_bo+flip(1-wage_domb_A).*QB_bo*(1+kapa)/((1+kapa)*(1+Par_b.trf_A))^(1/bet));

        qqbo_sqn = [sum(QA_bo)./QwA_bo; sum(QB_bo)./QwB_bo];

        wqbo_sqn = [(1-Par.bet)*Par.eta^(Par.bet-1)*qqbo_sqn(1,:).^(-Par.bet);
                    (1-Par.bet)*Par.eta^(Par.bet-1)*qqbo_sqn(2,:).^(-Par.bet)];  

        xkbo_sqn = [(Par.pi_frn/Par.bet)^(1/(1-Par.bet))./wqbo_sqn(1,:).^(1/Par.bet)*Par.L_B;
                    (Par.pi_frn_toUS/Par.bet)^(1/(1-Par.bet))./wqbo_sqn(2,:).^(1/Par.bet)*Par.L_A];

        II.IA_bo1_pr  = sum(prft_fncb_A.*QA_bo).*wqbo_sqn(1,:).^(1-1/bet);
        II.IA_bo2_rdi = sum(alfa_A/gama_A*Xb_iA.^gama_A.*QA_bo);
        II.IA_bo3_rde = sum(alfa_eA/gama_eA*Xb_eA.^gama_eA.*QA_bo);
        II.IA_bo4_wd  = pi_domb/(1-bet)*L_A*wqbo_sqn(1,:).^(1-1/bet).*sum(wage_domb_A.*QA_bo);
        II.IA_bo5_wf  = pi_frnb_toUS/(1-bet)*L_A*wqbo_sqn(2,:).^(1-1/bet).*sum(flip(1-wage_domb_A).*QB_bo);
        II.IA_bo6_wig = wqbo_sqn(1,:).*sum(QA_bo);
        II.IA_bo7_trf = Par_b.rev_US*Par_b.trf_A*((1+Par_b.trf_A)*(1+kapa)*eta/(1-bet))*wqbo_sqn(2,:).*xkbo_sqn(2,:).*sum(flip(1-wage_domb_A).*QB_bo);

        II.IB_bo1_pr  = sum(prft_fncb_B.*QB_bo).*wqbo_sqn(2,:).^(1-1/bet);
        II.IB_bo2_rdi = sum(alfa_B/gama_B*Xb_iB.^gama_B.*QB_bo);
        II.IB_bo3_rde = sum(alfa_eB/gama_eB*Xb_eB.^gama_eB.*QB_bo);
        II.IB_bo4_wd  = pi_domb/(1-bet)*L_B*wqbo_sqn(2,:).^(1-1/bet).*sum(wage_domb_B.*QB_bo);
        II.IB_bo5_wf  = pi_frnb/(1-bet)*L_B*wqbo_sqn(1,:).^(1-1/bet).*sum(flip(1-wage_domb_B).*QA_bo);
        II.IB_bo6_wig = wqbo_sqn(2,:).*sum(QB_bo);
        II.IB_bo7_trf = Par_b.rev_FN*Par_b.trf_B*((1+Par_b.trf_B)*(1+kapa)*eta/(1-bet))*wqbo_sqn(1,:).*xkbo_sqn(1,:).*sum(flip(1-wage_domb_B).*QA_bo);

        % Welfare
        % --------------------------------------------------------

        [WlfA,WlfB,IAb_set,IBb_set] = postestimation_welfare_trade(II,Par_b,horizons); 

        WLFA(:,ctr_s) = WlfA(:,1);
        WLFB(:,ctr_s) = WlfB(:,1);

        WLF_all(ctr_k,ctr_s).rateA = [trfA_call tau_bA_grid(ctr_s)];
        WLF_all(ctr_k,ctr_s).subsA = WlfA;
        WLF_all(ctr_k,ctr_s).subsB = WlfB;
        WLF_all(ctr_k,ctr_s).flag  = ind_neg_b;

    end    

    CEQ(:,ctr_k,:) = (WLFA./(Wlf_o*ones(1,size(WLFA,2)))).^(1/(1-psy));  

    WLFA_both(:,ctr_k,:) = WLFA;
    WLFB_both(:,ctr_k,:) = WLFB;

end

% save Optimal_rnd_over_open CEQ tau_bA_grid trfA_grid horizons 
% save Optimal_rnd_over_time_open CEQ WLF_all WLFA_both WLFB_both Wlf_o psy tau_bA_grid trfA_grid horizons McutX McutM MCUTX MCUTM

delete(gcp('nocreate'))
